Existence results for non-smooth second order differential 
inclusions, convergence result for a numerical scheme and 
application to the modelling of inelastic collisions 



o 

(N 



< 

U 



(N 
> 

o 
q 

o 
o 



X 



Frederic Bernicot 
CNRS - Universite Lille 1 
Laboratoire Paul Painleve 
59655 Villeneuve d'Ascq Cedex, France 
freder ic . bernicot @math. uni v-lille 1 . fr 



Aline Lefebvre-Lepot 
CNRS - Ecole Polytechnique 
CMAP 

91128 Palaiseau Cedex, France 
aline. lefebvre@polytechnique . edu 



March 9, 2010 



Abstract 

We are interested in existence results for second order differential inclusions, involving 
finite number of unilateral constraints in an abstract framework. These constraints are de- 
scribed by a set- valued operator, more precisely a proximal normal cone to a time-dependent 
set. In order to prove these existence results, we study an extension of the numerical scheme 
introduced in [5] and prove a convergence result for this scheme. 

Key-words: Second order differential inclusions ; Proximal normal cone ; Inelastic collisions ; 
Numerical scheme. 
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1 Introduction 

We consider second order differential inclusions, involving proximal normal cones. These ones 
were firstly treated by M. Schatzman [24] in the framework of elastic impacts and later by J.J. 
Moreau |13|. [T4"] to model inelastic impacts for a mechanical system in order to describe contact 
dynamics. The impact law describing the dynamics leads to a non- increasing kinetic energy at 
impacts. These second order problems appear in several models of mechanical systems with a 
finite number of degrees of freedom and dealing with frictionless and inelastic contacts. 

Let us specify this class of problems. Let I be a bounded time-interval, / : I x M. d — > M. d be 
a map and C : I M. d be a multi- valued map. The main question concerns the existence for 
solutions to the following second order differential inclusion: 

r Vi G I, q(t) G C{t) 
^+N(C(.),g(-)) 9 /(•,«(•)) 

Vt€J, q(t+) = P CtMt) q(t-) (1) 
?(0) =q Gint[C7(0)] 
, q(0) = uq. 
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We denote by int[C(0)] the interior of the set C(0), by N the proximal normal cone and for 
q 6 C(t), by Ct <q the set of feasible velocities: 

Ct,q '■= {u, q + eu G C(t + e) for small enough e > } . (2) 

We refer the reader to [3] and [lj for details concerning different normal cones ("limiting cone", 
"Clarke cone", ...). Here we will deal with "uniform prox-regular sets" C(t) so, according to 
[23J, all these cones coincide. 

Remark 1.1 We are looking for solutions q such that q has a bounded variation, in order that 
the second order differential equation in £7]) should be thought in the distributional sense. More 
precisely, we will solve it for time-measure q £ BV(I) and it should be written with time-measures 

dq + N(C(-),q(-))dtB f(;q(.))dt. 

In all this work, the second order differential inclusion will be written in the distributional sense 
for easiness. However, we emphasize that we consider time-measures. 

This differential inclusion can be thought as follows: the point q(t), submitted to the external 
force f(t,q(t)), has to live in the set C(t) and so to follow its time-evolution. The unilateral 
constraint ll q(t) £ C(t)" may lead to some discontinuities for the velocity q. For example, fric- 
tionless impacts can be modelled by a second order differential inclusion involving the proximal 
normal cone (see [131 14j). This differential inclusion does not uniquely define the evolution of 
the velocity during an impact. To complete the description, we impose the impact law 

Q(t + )=Pc tMt) q(n, 

introduced by J.J. Moreau in [13] and justified by L. Paoli and M. Schatzman in [T71IT9] (using 
a penalty method) for inelastic impacts. 

The set C{t) corresponds to a set of "admissible configurations" for q. In physical problems, it 
is generally described by several constaints (gi)i as follows : 

p 

C(t):=f]{q, gi (t,q)>0}. (3) 
i=l 

The existence of a solution for such second-order problems is still open in a general framework. 
The first positive results were obtained by M.P.D. Monteiro Marques [10] and L. Paoli and 
M. Schatzman [18] in the case of a smooth time-independent admissible set (which locally 
corresponds to the single constraint case p = 1 in ©). The proof relies on a numerical method 
involving a time-discretization of ([1]) in order to compute approximate solutions and is based on 
the study of its convergence . The multi-constraint case with analytical data was then treated by 
P. Ballard with a different method in [1], where a positive result of uniqueness for such problems 
was obtained. Then in [2U| ,l2H l2"2"]. an existence result is proved in the case of a non-smooth time- 
independent convex admissible set (given by multiple constraints). There, the active constraints 
are supposed to be linearly independent in the following sense: for each configuration q € dC, 
the gradients (V<7i(g))*gj associated to active constraints I := {i, gi{q) = 0} are supposed to be 
linearly independent. 

In the case of non-convex admissible sets, some results were obtained for a single constraint p = 1 
(for example in [B] or in [TT] and [2B] for the first result concerning time-dependent constraints). 
Recently in [8], B. Maury has proposed a numerical scheme for time-independent multiple and 
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convex constraints gi. The admissible set C is not supposed to be convex, however at each 
time step, the numerical scheme uses a local convex approximation of C. This improvement is 
interesting as it permits to define an implementable scheme, since the projection onto a convex 
set can be performed with efficient algorithms. 

A first result of convergence for this scheme was proved in [8] for a single constraint and appli- 
cations to the numerical simulation of sytems of particles submitted to inelastic collisions are 
studied in [7] by the second author. 

We emphasize that in the previously mentioned works, the different numerical schemes (per- 
mitting to discretize ([1])) are written (or can be written) in a multi-constraint case. The main 
difficulties consist in proving in a one hand the existence of solutions for (fTJ) and in the other 
hand a convergence result for the associated numerical schemes for such multi-constrained prob- 
lems. Concerning the uniqueness, we know from [53] and pQ that even with smooth data the 
uniqueness does not hold. The only positive results are proved in [25] for one-dimensional impact 
problems and in [1] in the context of analytic data. This critical question of uniqueness is not 
studied here. 

The framework 

In this work, we are interested in extending the previous work [5] , in order to prove the existence 
of solutions and to get a convergence result of the scheme in the case of multiple time-dependent 
constraints. Moreover we give applications in modelling inelastic collisions. 

First of all, let us precise some notations. We write W 1,oc (I, M. d ) (resp. W l,1 (I, W 1 )) for the 
Sobolev space of functions in L°°(/,R d ) (resp. L 1 (/, IR^)) whose derivative is also in L°°(/,lR d ) 
(resp. L 1 (/,R ,i )). BV(I,M. d ) is the space of functions in L°°(I, M. d ) with bounded variations on 
/. We define the dual space M(I) ■= (C C (J))' where C C (I) is the space of continuous functions 
with compact support (corresponding to the set of Radon measure due to Riesz Theorem). We 
set for the subset of positive measures. 

We consider second-order differential inclusions involving a set- valued map Q : [0, T] M. d 
satisfying that for every t G [0,T], Q(t) is the intersection of complements of smooth convex 
sets. Let us first specify the set-valued map Q. - This general framework has already been 
described by J. Venel in [25] for first order differential inclusions (fitting into the so-called 
sweeping process theory) and in [2] for a stochastic perturbation of such problems. - 
For i G {1, let g% : [0,T]xl rf ->Rbea convex function with respect to the second variable. 
For every t G [0, T], we introduce the sets Qi(t) defined by: 

Qi(t):={qeR d , 9i (t,q)> o} , (4) 

and the feasible set Q(t) (supposed to be nonempty) 

v 

Q(t):=C]Qi(t). (5) 

t=i 

We denote by I = [0, T] the time interval. The considered problem is the following one: we are 
looking for a solution q G W l >°°{I, R d ) , q G BV(I, R d ) such that 

' VtEl, q(t) G Q(t) 

|f + N(Q(-), <?(•)) 3 /(•,?(•)) 
' VtGJ, q(t+)=P CtMt) q(t-) ^ 

q(0) =q e int[Q(0)] 
, g(o) = n , 
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where N(Q(t),q(t)) is the proximal normal cone of Q(t) at q(t) and Ct >q is the set of admissible 
velocities : 

Ct, q :={u, d tgi (t,q) + (V q gi(t,q),u)>0 if gi (t,q) = 0}, (7) 
which corresponds to ([2]) in our framework. 

We have to make assumptions on the constraints g, t . First we require some regularity: we suppose 
that there exist c > and open sets Ui(t) D Qi(t) for all t in [0,T] verifying 

d H (Q i (t),R d \U i (t))>c, (AO) 

where dn denotes the Hausdorff distance. Moreover we assume that there exist constants 
a,/3,M > such that for all t in [0, T], gi(t, •) belongs to C 2 (Ui(t)) and satisfies 



and 



VqeUi{t), a< \V q gi(t,q)\ < P, 

Vq€Ui(t), \dtgi(t,q)\<p, 
VqeUi(t), \d t V qgi {t,q)\ < M, 
Vq€Ui{t), \B 2 q g t (t,q)\<M, 

Vq€Ui(t), \d 2 g t (t,q)\<M. 



(Al) 

(A2) 
(A3) 
(A4) 

(A5) 



In comparison with |28j and [2] where first order differential inclusion are studied, we require 
the new and natural assumption (|A5j) . due to the fact that we consider second order differential 
inclusions. 

Note that these assumptions can slightly be weakened. Indeed, the lower bound in (|A1[) is only 
required in a neighborhood of q G dQ(t). Moreover, we have assumed C 2 -smoothness in (|A4p 
and (|A5p for the sake of simplicity, but we only need C 1+e regularity. 

Furthermore, we require a kind of independence for the active gradients. For all t € [0, T] and 
q € Q(t), we denote by I(t, q) the active set at q 

I(t,q):={ie{l,..,p}, 9i (t,q)=0}, (8) 
corresponding to the active constraints. For every p > 0, we define the following set: 

I p (t,q):={ie{l,..,p}, gi(t,q)<p}. (9) 
We assume there exist 7 > and p > such that for all t £ [0, T], 



VqeQ(t), VAi >0, J2 X i\V q gi(t,q)\<7 

i€lp(t,q) 

We will use the following weaker assumption too: 



Vg € Q(t) , VAj > 0, *i|V g 0<(t,g)|<7 
iei(t,q) 



Xi^ q 9i(t,q) 



XiV q gi{t,< 

iel(t,q) 



(A6) 



(A6' 



Note that assumptions (|A6p and (|A6'P describe a kind of "positive linear independence" of the 
almost active gradients. In the time-independent case, (|A6'P is lightly weaker than the linear 
independence assumption made in [20\ |2~T| |2"2] , In fact, such assumptions imply a "uniform 
prox-regularity" of the admissible set Q(t), which is a weaker property than the convexity. 

Under these assumptions, we have a characterization of the proximal normal cone N(Q(t), •). 
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Proposition 1.2 (Prop 2.8 of |28j) In this framework, we know that for every t G I, and 
every q G dQ(t), 

N(Q(t),g) := < - x * ^ \ ■ 

{ iel{t,q) J 

So our Problem ([U]) can be written as follows: we are looking for solutions q G W l, °°(I, 
BV(I,M. d ) and time-measures Aj G .M+(I) such that 

' Vt G J, g(i) G Q(t) 

2, P 



/(•,?o) + E^ft(-«(')) 



i=l 



supp(Aj) C {t , gi(t, q(t)) = 0} for all i 



(10) 



Vt G I, q(t~ 



c. 



t,q(t) 



q(0) =q e int[Q(0)] 
. 9(0) = uq. 



We denote by A = (Ai, 
p constraints. 



X p ) G W the vector of the Lagrange multipliers associated to these 



As usual, we obtain existence results for (|1Q[) by proving the convergence of a sequence of dis- 
cretized solutions. 

To do so, we extend the algorithm proposed by B. Maury in [8] for modelling inelastic collisions, 
to the case of abstract and time-dependent constraints. In [8j, the convergence (up to a subse- 
quence) is proved in the case of a single constraint. Here, we show that this convergence still 
holds in the multi-constraint case. 

Let us describe the numerical scheme. Let h = T/N be the time step. We denote by q^ G M. d 

and u 1 ^ G M. d the approximated solution and velocity at time Viji for n G {0, ..,iV}. 

The discretization of the continuous constraints Q« )g (t«) proposed in [8] corresponds to a first 
order approximation of the constraints in the space variable: for t 6 7 and q G U(t), we set 



Kh(t,q) := {u, gi(t,q) + h(V q gi(t,q),u) > 0} 



(11) 



The reason why we do not expand the time variable in this discrete admissible set is that we 
will use a semi-implicit numerical scheme, with directly impliciting in time the set Kh(t,q) (See 

(USD). 



The approximated solutions are built using the following scheme: 
1. Initialization : 

(<lh, u h) : = («D,«o) 



(12) 



2. Time iterations: q 1 ^ and are given. We define := — 



f(s,q%)ds, 



u 



n+1 . 



and 



:=ql + hul+\ 



(13) 
(14) 



5 



where Pc is the Euclidean projection onto the set C. This algorithm is a "prediction-correction 
algorithm": the predicted velocity + hffi, that may not be admissible, is projected onto the 
approximate set of admissible velocities. 

Since the projection ,. n +i n \ consists in a constrained minimization problem, with a finite 

number of affine constraints, it involves Lagrange multipliers (A^ +1 ) G W corresponding to the 
p constraints. It can be checked that we have a discrete counterpart of the momentum balance 
appearing in (fT0|) : 

= fk+E Kf^Mtl + \€) (i5) 

with \ n h f > and A"+ x = when g t (t n h +1 ,q^ + /i(V ? > 0. 

In [8], the scheme is shown to be stable, robust and to present a good behaviour for large time- 
steps. That is why, we are interested in continuing its numerical analysis in the multi-constraint 
case, with proposing some extensions like the time-dependence of the constraints. 

Results 

We recall that / = [0, T] is the time interval and h is the constant time step ifth for n = . . . N). 
We denote by the piecewise affine function with qhifih) = Qh- We denote by Uh the derivative 
of qh, piecewise constant equal to u^ +1 on Finally, we define Xh, piecewise constant 

equal to A™ +1 on ]ijj,tj| +1 [. 

The convergence theorem is the following one. 

Theorem 1.3 Let (qhi u hi^h) be the sequence of solutions constructed from the scheme Iil2]\15\) 
and suppose that f : I xM. d — >• M. d is a measurable map satisfying: 

3K L >0, VtG J, Vq,qeU(t), \f(t,q)-f(t,q)\<K L \q-q\ (16) 
3F e L\I) , Vt G /, Vg G ?7(t) , |/(t,?)| < (17) 
Then, when h goes to zero, there exist subsequences, still denoted by (qh)h, (uh)h, (^h)h, and 

(q,u,X) G W l '°°(I,M d ) x BV(I,R d ) x M+{I) P 

such that 

u h — >u in L 1 (/,R d ), 

q h — >q in W 1 ' 1 ^,^) and L°°(/,]R d ) with q = u, 

X h ^Xin M+{I) P 
where (q, u, A) is solution to [W\) and so (q, u) is a solution to 

We emphasize that, up to our knowledge, this result is the first one concerning such multi- 
constrained second order differential inclusions with on the one hand time-dependent constraints 
and on the other hand a non-convex and non smooth admissible set. 

Remark 1.4 For time-independent constraints, Assumption iA6)) is required but Assumption 
HA6\) is not necessary. 

The proof is quite long and technical. We refer the reader to [8] for a first proof dealing with 
the case of one (p = 1) time-independent constraint g. We will follow the same reasoning with 
some new arguments (appearing in [28] ) in order to solve the difficulties raised by the multiple 
constraints and the time-dependence. Section [2] is devoted to the outline and the main ideas 
of the proof. For the sake of readibility, the demonstrations of some technical Propositions are 
postponed to Section [3J In Section HI we describe an application to the modelling of inelastic 
collisions. 
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2 Convergence result 



This section is devoted to the proof of Theorem 11.31 It is divided in 7 steps and for readibility 
reasons we have postponed some technical proofs in the next section. 

• Step 1: The scheme is well-defined and produces feasible configurations 

Proposition 2.1 For a small enough parameter h, the scheme is well-defined. Moreover the 
computed configurations are feasible : 

V/i>0,VnG{0,..,iV}, qh {tl)eQ(tl). 

Proof : Let h be smaller than c/cq where c and cq are given in Lemma 12.21 (below stated) 
and Assumption (jAOp respectively. By assuming that qhit^) G Qifih), we also deduce that 
q h {tl) e Q{tl +l ) + c hB(0, 1) C U{tl +l ). Then the gradient V q gi (tl +1 , q%) is well-defined and 
so is the set , gJJ). The step 2 of the scheme can be performed and due to the convexity 

of function 0> 

ul+ l eK h {tl+\ql)^ql+ l eQ{tl +l ). 
Then we conclude by iteration. □ 

Lemma 2.2 The set-valued map Q is Lipschitz continuous with a constant cq, for the Hausdorff 
distance. 

We refer the reader to Proposition 2.11 of [28J for a detailed proof of this result. 

For the intermediate times t the point qh(t) may not belong to Q(t). However from 

Proposition 12.11 and Lemma |2.2| we have the following estimate : 

Vfc>0,Vi€l, d(q h (t),Q(t)) <^{d{q n h ,Q{t)),d{q n h + \Q{t))} < c h. (18) 

• Step 2: u h is bounded in BV(I, R d ) 

First, we check that the velocities are uniformly bounded (proved later in Subsection 13. ip . 
Proposition 2.3 The sequence of computed velocities (uh)h is bounded in L°°(I,M. d ). We set 

K := sup 1 1 Uh < 

h 

Proposition 2.4 The sequence of computed velocities {uh)h is bounded in BV(I,M. d ). 

Proof : Since u° h = uq and using Proposition 12. 3\ it suffices to show that (uh)h has bounded 
variations on I. This has been proved for a single constraint in [8J. Unfortunately, this proof 
cannot be extended to the multi-constraint case. To obtain an estimate on the total variation, 
we use a similar technique to the one proposed in [5j and [6j. These ideas rest on the following 
property: all the cones K^t^ 1 , qfi) contain a ball of fixed radius with a bounded center, which 
describes the fact that the solid angles of the cones N(Q(£), qh(t)) are not too small. This 
property is proved by using a "good direction" (see Lemma 13. ip which permits to increase all 
the almost active constraints. The details of the proof are postponed to Subsection 13.11 □ 

• Step 3: Extraction of convergent subsequences 

Proposition 12.41 directly implies the following convergence result: 
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Proposition 2.5 There exist q in W 1)OC (I, M d ) and u in BV(I,R d ) such that, up to a subse- 
quence, 

u h ► u in L 1 (/,R d ), 

h— >-0 

% ► q in W 1 ' 1 (J> K<1 ) and L°°(I,R d ) with q = u. 

Furthermore U8\) yields 

Vt € /, g(t) 6 Q(i). 
In addition, we show that the sequence of Lagrange multipliers converges: 
Proposition 2.6 There exists A in A4+(I) P such that, up to a subsequence, 

X h A in M+{I) P . 

Proof : From (|15|) . we have 

h\ n h fV q 9l {tl + \ql) = < +1 - ul - hfl 

i 

with A™ +1 > and A™ +1 = when 

gi (tl + \ql) + h(V qgi (tl + \q n h ),ul +1 ) >0. (19) 

Remark that for a small enough parameter gi(t^ +1 , q%) + h{V q gi(t^ +1 ,q^),u]^ + /i/^) < 
implies that i E 2p(£^ +1 , g^). Consequently, the reverse triangle inequality (Assumption (|A6'|) ) 
with (|Aip and the Lipschitz regularity (Assumption (|A3P ) imply for a small enough parameter 
h 

h X n+1 < \n n+1 n n hf n \ 

11 h,i ~\ U h U h n Jh\- 

i 

Therefore, we obtain for this small enough parameter h 

UhWmi) ^ Var 7 (u fe ) + IIFHii^). 

Hence from Proposition 12.41 together with hypothesis (|17p . (Xh)h is bounded in L l (I), which 
concludes the proof. □ 

• Step 4: Momentum balance 

As in Step 6 of Theorem 1 in [8] , using Proposition 12.51 together with Proposition 12.61 we pass 
to the limit in the discrete momentum balance (1151) to obtain 



Proposition 2.7 The momentum balance is verified by the limits u and A: in the sense of 
time-measure 

U = /(■,<?(•)) + X iV q gi(;q(-)). 

i 

Note that this equation has to be thought in term of time-measure (see Remark : 

du = f(-)dt + Y^V q 9i(;q(-))k, 

i 

where we denote by du the differential measure of the l?l/-function u. 
• Step 5: Support of the measures Aj 

From the uniform convergence of q^ and the Lipschitz regularity of gi (Assumptions (jAip and 
2j)), it can be checked, as in Step 7 of Theorem 1 in [8], that 
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Proposition 2.8 

Mi , supp(Xi) C {t, gi(t,q(t)) = 0}. 

Indeed (|19|) describes a similar property for the discretized multipliers. The uniform convergence 
allows us to go to the limit in (fT9|) and to prove the previous proposition. 

This property describes the fact that the measure Aj has a contribution only when the associated 
constraint gi is saturated. 

• Step 6: Initial condition 

As in Step 8 of Theorem 1 in [8j, using again the uniform convergence of qh it can be shown 
that 

q(0) = qo and ti(0) = uq. 

We emphasize that to prove this point, we use the property go £ Int[Q(0)]. From this, it can 
be checked that for < s with s a small enough parameter the desired velocity u% + hf'^ still 
remains admissible and so we don't need to project. That allows us to deal with any initial 
velocity uq G If qo G dQ(0), this property still holds if we assume that the initial velocity is 
admissible: uq G Co, go . Else, we would get 

^ + (0) = Pe„, 9 >o) 

according to the next proposition. 

• Step 7: Collision law 

Finally, Theorem II .31 will follow, provided that we check the collision law for the limits u and q, 
which is given by the following proposition. 

Proposition 2.9 

VtoGl, u + (t ) = F CtQMtQ) (u-(t )). 
Proof : The idea is to let h go to zero in the discrete collision law, 

The main difficulty comes from the fact that the mapping q — > Kh(t, q) is not Lipschitzean. The 
details of the proof are postponed to Subsection 13.21 □ 

3 Auxiliary results 

Before proving the technical Propositions 12.31 12.41 and 12. [)\ we recall the following main Lemma. 
This technical result is very important and all our proofs rest on this idea. It says that, for each 
time one can find a "good direction" increasing all the constraints which are almost active, 
the corresponding increase being independent of n. 

Lemma 3.1 There exist constants 5,k,0 and r > such that for all t £ I, for all q € dQ(t) 
there exists a unit vector v := v(t, q) satisfying : 

• for all s G [t — r, t + t], y G Q(s) n B(q, 6) and i G I K p(s, y), 

(V q gi{s,y),v) > 5, 

where p is the constant defined in iA6\) . 

This lemma is a consequence of the reverse triangle inequality (Assumption (|A6|) ). We do not 
give the proof here and refer the reader to Lemma 2.10 in [28] for a detailed proof with r = and 
to Proposition 4.2 in [2J for a complete proof with some r > 0. Indeed there is in the previously 
mentioned papers, a detailed construction of such "good directions". 
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3.1 BV estimate for u h (Propositions [2731 and I2T41) 

We first prove a uniform bound of the computed velocities Uh in L°°(I). 
Proof of Proposition 12.31 

1 — ) For any t G I and g G Q(t), construction of a specific point w G Ct i(? . 

From Lemma I37L1 (stated at Section [3|), it exists a "good direction": a unit vector v satisfying, 

Vielfaq), (V qgi {t,q),v)>8, (20) 

for some numerical constants S, 1 > non-depending on t and q. For k > a large enough real, 
we claim that few belongs to Ct. q . Indeed for i G Ii(t, q) (d /(t, g)), with fc > (/3 + S)/6, it comes 

dt9i(t,q) + (V q gi(t,q),kv) > d t9i (t, q) + kS > S > (21) 

thanks to Assumption (|A2|) . 

Choosing k = (ft + 5)/5, we have built a point w = kv belonging to Ct s with |w| = (f3 + 5)/5. 

2— ) This vector u> belongs to Kh{s + /i, q) for (s, g) close to (t, q) and /i small enough. 

Let fix a point (t, q) (for example the initial condition (t, q) = (0, go))- From the previous point, 
we know that there exists a bounded admissible velocity w G Ct t q, which satisfies the stronger 
property ([2Tjh 

More precisely if a € I and q G Q(s) satisfy 

\s-t\ + \q-q\< min{//(2/3), e} := v (22) 

with e a small parameter verifying 2eM (1 -\-(ft + 5)/5) < (5 then u; G i£/i(s + ^ 5 ?)> Indeed, from 
(1221 we have 

h/2{s,q) C IK*)?)- 
This, together with ([2T]) and ([22]) gives, for all i G J^/ 2 (s, g) 

+ <Vg0i(s, £),«;) >*-eAf(l + H) ><5/2. 

Moreover, for the other indices i ^ h/2{ s -,Q) we have 

&(s,g) +/i[%i(s,g) + (Vgfft(s,g),w)] > ~ - ^/3(1 + 
Consequently, for /i small enough (/i < /to := 1/(2/3(1 + (/3 + 5)/<5) + 5/2)), we obtain 

Vi, 5i(s,g) + h[d t gi(s,q) + (V g t/i(s,g),«j)] > /i-, 
which, by a first order expansion in time gives: 

Vi, gi(s + h,q) + h{\7 q gi(s + h,q),w) >h- + O h ^ (h 2 ). 

We deduce that there exists h\ < ho such that, for h <h\, 

Vi, gi(s + h,q) + h(V q g^s + h, q),w) > 0, 

and consequently w G K^s + h,q). 

3— ) Estimate on the velocities for small time intervals. 

Let us fix h < hi (given in the previous point) and a small time interval [t_,t_|_] C I of length 

It i < V 

U ~ t ~ ~ 2(\ul°\ + 2^ + f T F(t)dt) 
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where no is the smallest integer n such that > i_. We suppose that G [t_,i+]. We are 
looking for a bound on the velocity on this time interval. From the first two points, setting 
(t,q) = (^°,<Z^ ), we have an admissible velocity w G C^ q such that for all s G I and q G Q(s) 
we have 

w G K/i(s + h,q) 

as soon as |s — 1\ + |g — q\ < za Since for s = G t + ] , \s — t\ < v/2, we deduce that for all 
integer n such that t% G [£_,£+], if 

\q n h -q n h °\<u/2 (23) 



weK h {tl + \qt) 



then 

Considering such an integer n satisfying (f2"3"j) . since is the Euclidean projection of u^ + hfj^ 
on the convex set Kh(t^ +1 , g^) (containing the point to), we deduce that 

\ul+ 1 -w\<\ul + hfZ-w\ 1 

which implies 

- w\ < \ul -w\ + / F(t)dt. 



We set m, the smallest integer (bigger than no) such that m + 1 does not satisfy (f2"3"j) or t™ +1 ^ 
[£_, £+]. By summing these inequalities from no to n = p — 1 with no < p < m, we get 

VpG[n ,m], |it£ < + /" F(£)d£. 

•/ii 

Finally, it comes 

sup \ u p h \<\ul°\ + 2^-+ [ F(t)dt. (24) 

no<p<m JO 

By integrating in time, we deduce 

- CI < (k°I + 2 ^ + £ F ^ dt ) i*+ - *-i ^ 

by the assumption on the length of the time interval. As a consequence, we get that n = m + 1 
satisfies (j23|) which by definition of m, yields < t+ < tV? . Hence, from (|24|) . we have 



sup |uJJ| < |u"°| + 2^±i + f T F{t)dt. 

t-<ti<t+ ' o Jo 



4—) End of the proof. 

The parameter h < hi being fixed, we are now looking for a bound on on the whole time 
interval I = [0, T]. Let us start with £_ = £(0) := 0. From the previous point we know that with 

t+ = t(l) := min < 

we have 



2(\u \ + 2^ + J T F(t)dt 

sup \u h {tt)\ < \u Q \ + 2^-^ + [ F(t)dt 
t"<t(l) o Jo 



o<q<t(i) 
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Then, let us suppose that there exists m such that t(0) < t 1 ^ < t(l) < t^ 1+1 . We have 
< Si := t(l) - t^ 1 < h. In that case, we set i_ = t^ 1 and 



t + = t(2) := min { t ni + 



mm < 



2(\u \+^ + 2 jT F{t)dt 
- 5i + 



,T 



2(\u \ + A^ + 2^F(t)dt 



From the previous point, we deduce that 



sup |u h (^)|< sup \u h (tl)\ < + 2 / F(t)cfc 

i(l)<t£<t(2) <" 1 <t?<<(2) " jO 



t£<i(2) 



and so 



sup \u h (tl)\ < \u \ + + 2 / F(i)di. 

:t"<t(2) Jo 



0<*™<t(2) 

By iterating this reasoning, for any integer k > 1 we set 
t(k) := min <J t(k - 1) - + 



2 nw, | + 2A:^ + F(t)dt) 



mm 



> . 



i=i .=1 2 [\u \ + 2i^ + 2i F{t)dt) 
where Sk < h for all k. This construction of t(k) can be made while there exists nj~ such that 



t(k - 2) < tl"- 1 < t(k - 1) < That is, while i(Jfc - 1) - t(jfc - 2) > /i. This condition will 



be verified as long as 



"4-2 + 



2 ( |uq| + 2(fc - + 2(k - 1) Jq F(t)dt 



T 



Therefore, using the fact that < < h, we see that we can construct t(k) for k < N 

verifying 

;/ 

. > 2/i, 



(25) 



2 M«o| + 2(k - l)2+£ + 2(jfe - 1) F(t)dt 
which is equivalent to 

*<^^-¥)(^+fH _1 - 

Consequently, we know that the velocities can be bounded on [0, t(ko(h))] where 



ko(h)-l ko(h) 

t(ko(h)) = min I- V <5i + V — t — • 

4=1 i=i 2 K| +2i^±i + 2i/ / F(t)di 



Now, using the fact that fco(/i) goes to infinity when h goes to zero, that the harmonic serie 
diverges and that ([25]) yields 



fc-1 

i=i 



< /ifco(/i) < C, 
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we see that t(ko(h)) is equal to T for h small enough. Therefore, there exists h 2 < h\ such that, 
for h < hi-, T = t(ko(h2)) = t(ko(h)). Finally, we see that, for h < h 2 , t(k) can be constructed 
until k = ko(h 2 ) and Uh can be bounded as follows 

sup sup \u h (tl)\ <\u \ + 2ko(h 2 )^^- + 2k (h 2 ) [ F(t)dt (26) 

'i<h 2 0<t"<T Jo 



h<.}%2 u^c^ 

which concludes the proof of the existence of a uniform bound in L°° for the velocities u^. □ 
To prove Proposition 12. 4| it suffices now to show that the sequence (v,h)h has bounded variation. 

Theorem 3.2 The sequence (uh)h has a bounded variation on I. 

Proof : In order to study the variation of Uh on I, we split I into smallest intervals. We define 
(sj)j for j from to P such that: 

s = , s P = T 

\Sj+l 



T, 




1 . / 




— mm < 
2 I 


4} 


1 




- — mm 
" 2 


B 



where r and 6 are given by Lemma 13. II and K is the bound on (see Proposition 12 .3[) . 

All these constants do not depend on h and such a construction gives 



P 



2T 



mm 



+ 1, 



(27) 



which is independent of h. Then, for all h, we define n 3 h for j from to P — 1 as the first time 
step strictly greater than Sj\ 

and nf{ is set equal to N (t^ = t h h = T). 

In the following, we suppose h < min{|sj + i — Sj\}/2. Doing so, we obtain a strictly increasing 

j 

sequence of (t h h )j with 

(28) 



|t. " — r ft I < mm < r, 

The variation of on I can be written as follows 

AT-l 

Var/(u fc ) = £ K +1 - ul\ = 

n=0 



A' 



p-i 

i=o 



where 



n{ +1 -l 



Var,K):= £ 



corresponds to the variation on [t^ h , [. To study these terms, we recall that 

< +1 = P K h{ t^, qV K + hf n ] 
by construction and state the following lemma: 



(29) 
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Lemma 3.3 There exist rj > and uniformly bounded vectors y n h such that, for all small 

j j4 



enough h, for all j = . . . P and n £ [n{ , n{ +1 [ , to e have 



-n r n i I 1 f i \2 \ i2 

si = P x h (^ +1 ,^)NJ =>• \%i - x \ < — [\x - y h\ -\ Xl -yh\ 
Proof : The outline of the proof is the following: first, we prove that there exist unit vectors 



v n h such that 

n £ [nin{ +1 [ =► B(^fv<, V ) C K h (t n h + \q n h ) with r, : = |, (30) 

where if is a bound on Hn/JI^^) (see Proposition 12. 3p . Then, we conclude using similar 
arguments to the ones exposed in [6] . 

j 3 

Step 1: From Lemma 13.11 with t = t, h and q = q h h , we have a unit "good direction" written 
v n h. Let n belong to [n^,n^ +1 [. From Proposition 12. 1| we know that <7^ +1 belongs to Q(^ +1 ). 

7 7 

Moreover, ([28]) gives — < r and — q h h \ < 6. Consequently, Lemma 1331 gives 

Vi e/X 1 .^ 1 ). (Vgft^.g^.X) >& (31) 
We deduce that for all index i £ {1, ..,p} and a small enough parameter h : 

ftC^.^ + ^MV^C^.g^ 1 ),^) > -hf3K + 2h(3K = hf3K. (32) 
Indeed, we write 

9i(tt\ql) + ^h(V qgi (t n h + \q n h +1 ),v<) = 

{ 9l (t n h + \q n h ) - 9i(t n h + \q n h +1 )} + ^(t n h + \q n h +1 ) + ^fh{V q9i (tl+\q^ 1 ),V<^ . 

The first term can be estimated using (jAip and the bound if on ||'Uh||.L°°(n- In order to estimate 
the second term, if i £ J Kp (i£ +1 , q% +1 ), we use (JHTJ) together with the fact that gi(t% +1 , q% +1 ) > 0, 
which gives the required bound. In the case i £ J re p(t? + , g^ +1 ), we use ) > ftp 

and (jAip . which also gives the required bound for /i small enough. 
Finally, ([32]) together with assumption ([AT]) implies that for all u £ B(^^v n h,K/2) 

9l (t n h + \q n h ) + h{V qgi (tl + \q n h lv) > h ^K-^ = ^f>0, 
which proves (|30j) . 

Step 2: Let n belong to [n^,n^ +1 [. We define 

z n h := y n h + r/— i- where y n fc := — — t/V 

F0 — xi\ o 

(Here we suppose xo ^ x\, else the desired result is obvious.) From the previous step we have 

AeB(^v<,r,)cK h (tl +l ,q%). 
The point x\ being the projection of xq onto the closed convex set K^tV^ 1 , q 1 ^) , we get 

(x - Xl, Z n h - Xl) < 0. 
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From this we have 

\x - y nJh \ 2 = \xi - y<\ 2 + \xq - xi| 2 + 2{z< - y<,x - Xl ) + 2{ Xl - z<,x - x x ) 

> \x l -y<\ 2 + 2(z<-y<,x -x l ) 

> \xi - y n h\ 2 + 2r]\x - x{\. 

j q a jy~ 

This, together with the fact that the vectors y n h are uniformly bounded by -t— , ends the proof 
of Lemma 13,31 □ 

We now come back to the proof of Theorem l3.2l For n in [ni, using (|29|) and the previous 

lemma (with xq = u\ + hffi and x\ = ), it comes 



lu^-ul-hfZl < i^o-y^l 2 -^-^! 2 ) 

1 /. - - - n a n{\2 L.n+1 |2 



where L := 2/3K/5 (see Proposition 12.31 for the definition of K). By summing up these terms 



for n from ni to n 1 ?" 1 " 1 — 1 we get 



2 



Var,K)= £ K +1 -<l < ^(K 1 -^i 2 -K r -^| 2 )+ £ 

+ I ( ^ + L + r? ) £ 1/,/nj 



and finally 

P-l , P-i 



1 / 3 3+1 \ 1 

varK) = ^v arj M < ±Y,(\< h -y ni \ 2 -K h -y ni \ 2 ) +^W F W 



2 

Ll(I) 



3=0 ' 3=0 

1 



+-(K + L + r,)\\F\\ L1{I) 
< I ( K + L) 2 P + + J (K + L + r,) \\F\\ LHI) . 

This completes the proof of Theorem 13.21 since P does not depend on h from ()27p . □ 
3.2 Collision law for the limits u and q (Proposition 12. 9|) 

This subsection is devoted to the proof of Proposition 12. 9\ recalled in the following Theorem : 
Theorem 3.4 Let to £ I be fixed. The limit function u verifies: 

u+(to) = P Cto , q(to) (u-(t )). 
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Note that, from Proposition 12.51 u G BV(I), so that the left-sided u~(io) and the right-sided 
u + (to) limits are well-defined. 

The proof is quite technical so for an easy reference, we remember the definitions of the sets Ct, q 
(given in ([?])): 

Ct,q ■■= {u, d t gi(t,q) + (V g gi(t,q),u) > 0, if <?) = 0} 
and K h (t,q) (given in (QT])): 

K h (t,q) := {u, gi(t,q) + h{V q 9i(t,q),u) > 0} . 

Moreover, we recall that 

K := sup 1 1 ith ||£oo(7) < oo. 

h 

The desired property 

u + (t ) = P Ctoq(t Ju-(t )) (33) 
can be seen as the limit (for h going to 0) of the "discretized property" 

< +1 = ?K h (tl+\ q vl< + hn (34) 

Proof : First we claim that 

u + (t ) G C t0A[tQ) . (35) 

To verify this property, let us consider an index i such that gi(to,q(to)) = 0. Then a first order 
expansion gives : 

9i(h + e,q(t + e)) = e [d t gi(t , q{t )) + (u + (t ), V q gi(t , q(t )))] + o e _> (e)- 
The feasibility of q(to + e) (see Proposition 12. 5p yields 

dt9i(to,q(to)) + (u + (t ), V qgi (t ,q(t ))) > 
which corresponds to ([35]) . 

Let us now come back to the proof of (f33|) . As we just proved ii + (to) G Q ,g(*o) anc ^ smce Q ,g(*o) 
is a convex set, ([33]) is equivalent to 

Vu; GC Mto) , (n^(to)-n + (t ),^-w + (t )) < 0. (36) 

So, in the following, let us choose w G C io )9 (t )- To prove (f36|) . we construct a family of points 
w v for ^ > such that w v tends to w when v goes to zero and satisfies w v G Kh(t + /i, q) for h 
sufficiently small and (t, q) close to (to, q(to)). Then, for each v, we go to the limit on h, t and q 
to show that (u~(to) — u + (to),w u — u + (to)) < and finally, we make u go to zero to conclude. 

Step 1: From Lemma [3.11 there exists a neighborhood U C I x~R d around (to, q(to)) and v € M rf 
such that for all t G / and q G Q(t) 

(t,q)eU VieI Kp (t,q), (V, gi(t,q),v) >5, (37) 

with a numerical constant J > 0. For ^ > 0, we consider the point to„ := io + z^u with > 0. 
For all i G I Kp (t,q) H I(t ,q(t )), dSTJ) together with to G C to , g (t ) gives 

dtgi(t ,q(t )) + (V q gi(t,q),w v ) = d t gi(t , q(t )) + (V q gi(t,q),w) + u(V q gi(t,q),v) 

> (Vq9i(t,q) - V q gi(t ,q(t )),w) + za5 
>i/*-AfH[|t-to| + |9-9(to)|] 
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and consequently from Assumptions (jA3J) and (|A5p 

d t g i (t,q) + {V q g i (t,q),w„) > vS - (M\w\ + M) [\t - t \ + \q — q(to)\] ■ 

So for every v > 0, if (t, q) is closed enough to (io><z(*o))> we deduce that for all i S I Kp (t,q) n 
I{t Q ,q{t Q )) 

9t9i(t,q) + {V q gi{t,q),w v ) > —. 
For the indices i ^ <?), we have 

5i(t,g) + h[d t gi(t,q) + (V q gi(t, q), w v )\ > np - hf3(l + + i>). 
Finally for j ^ I(t ,q(t )), 

9i {t, q)+h [d t gi(t, q) + (V, gi(t, q),w u )] > a - h/3{l + |w| + v) - [\t -t \ + \q- q(t )\] , 
with 

a := min gi(t ,q(t )) > 0. 
i$I{to,q(to)) 

We conclude that for each fixed v > 0, there are and /ij, such that for every h < h v , (t,q) € U 
with \t — to I + \q — q(to)\ < e ^ an d 9 € Q(t) we have 

Mi, gi(t,q) + h[d t gi(t,q) + (V ? gi(t,q),w v )] > h—, 



which by a first order expansion in time gives: 



vS 



Vi, + /&,<?) + h(V q gi(t + h,q),w v ) > h— + O h ^ (h). 

At the cost of decreasing h u , it comes for h < h u , 

Vi, gi(t + h,q) + h(V q gi(t + h,q),w v ) > 0, 

and consequently, wv 6 i£ft(t + /i, g) for every h < h v , (t, q) € U with \t — to \ + \q — q(to)\ < e u 
and q E Q(t). 



Step 2: Let us now fix the parameter v. 

Thanks to the uniform Lipschitz regularity of the maps q^ and their uniform convergence towards 
q, there exists h v < h u such that for e < e u /(2 + 2K) and h < h u , 



t£,t£ +i € [t - e,t + e] => [t£ +i -t \ + \q k h - q(t )\ < e v . 

From this, as q\ S Q^X) from Proposition 12. 1\ the previous step (with t = t|) gives G 
Therefore, Kf l (t^ l +1 , q%) being convex, we have 

(4 + hf k h -u k h + \w„-u k h +1 )<0. (38) 

We sum up these inequalities for k from n to p, integers chosen such that tV: is the first time 
step in [to — e,to — e + h] and t^ +1 the last one in [to + e — h, to + e]. First, we know that 



j2h{f\w u -u k+1 ) 



< (Kl+if) [° + F(t)dt 

J tn— e 



(39) 
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with K := supfj ||uh||oo- We also have 



J>£ - u k h + \w u ) = {u h (t n h ) - u h (tl +1 ), Wv ). 



(40) 



We deal with the remainder as follows: we write 

p p 



«|2 i |„P+1|2 
ft I ' 



which gives 



k+l\2 



l r 



+ 2 



P+l\|2 



^Var 2 («fc)L,tf] + | [-|^(^)| 2 + K(*DI 2 ] . 



(41) 



where we wrote Var 2 for the L 2 -variation of a function. Using ([38]) . ([39]) . (|40|) and (|4T|) . we 
finally get : 



\\^ 2 {u h f KK] + ~ [-\u h (tl)\ 2 + \u h (tl +1 )\ 2 ] + {u h (tl)-u h (tl),w v ) < [°^F(t)dt- 



Let us now choose a sequence of e m going to zero, such that Uh pointwisely converges to u at the 
instants to — e m and Iq + e m (which is possible as Uh converges almost everywhere towards u) . 
For each e m and h small enough, we have shown that the last inequality holds. Then, passing 
to the limit for h —> we get 



9 Vax 2(«)ft _ Wo+ero ] + \ [-\u(t - e m )\ 2 + |t*(*o + e r , 



+ (u(t Q - e m ) - u(t + e m ),w v ) 



< 



t()—£r, 



F(t)dt, 



which gives for e m — > 

^Var 2 {u) 2 t - )t + ] + -[-\u-(to)\ 2 + \u + (to)\ 2 ] + («~(i ) - u + (t ),w u ) < 0. 
Finally we obtain 

- \u + (t ) - u-{t )\ 2 + - [-\u-(t )\ 2 + \u + (t )\ 2 ] + (u-(to) - u + {t ),w u ) < 0. 
By expanding the square quantities, this can be written as follows 

(u-(to)-u + (to),w v -u + (t )) < 0. 



(42) 



To conclude the proof, it now suffices to remember that w v = w + vv and, since for each v > 0, 
the previous reasoning holds, we obtain by letting v go to in ([42 j) . □ 
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Figure 1: Particles i and j : notations. 



4 Application to the modelling of inelastic collisions 

The continuous model 

We consider a mechanical system of N spherical rigid particles in three-dimensions. We denote 
by qi G M 3 the position of the center of particle i, by ri its radius, by mi its mass and by /, £ M 3 
the external force exerted on it. Let q 6 M. 3N be defined by q :=(..., g^, . . .) and / E R. 3N by 
/ := (. . . , fi, . . .). We denote by Dij(q) the signed distance between particles i and j: 

Ai(g) := h -qj\- in + rj), 
and we set eij(q) = (qj - qi)/\qj - qi\ (see Fig. [TJ. 

The problem we are interested in is to describe the path of the configuration q submitted to the 
force-field / and undergoing inelastic collisions. This inelastic collision law can be modelled by 
imposing non-overlapping contraints on the particles (see the work of J.J Moreau [15J introducing 
this concept). Therefore, we write that the positions of the particles have to belong to a set of 
admissible configurations Qo avoiding overlappings: 

?gqo : =n^' Ai(<?)>o}. 

id 

We define M as the mass matrix of dimension 3Nx3N, M = diag{. . . ,mi,mi,mi, . . .). Then, we 
denote by Gij G M. 3N the gradient of distance Dij with respect to the positions of the particles: 

Gij{q)= (...,0, -dj ,0, ...,0, aj ,0,...,0)*. 

* 3 

The set C q is the set of admissible velocities: 

C q := {u, (G i:j (q),u) > 0, if Dijiq) = 0} . (43) 

To finish with notations, we denote by A = (. . . , Ay, . . .) € M^^ -1 )/ 2 the vector made of the 
Lagrange multipliers associated to the N(N — l)/2 constraints u Dij(q) > 0". 

Let / =]0, T[ be the time interval. The multi-particle model we are interested in may be formally 
phrased as follows: 

' q € W 1,00 (I,R m ) , q € BV{I,R 3N ), A € {M + {I)) N{N - 1)/2 , 

Vtel, q(t + )=P Cqit) q(t-) 

Mq = f + J2\ ij G ij (q) 

i<j (44) 
supp(Ay) C {t , Dij(q(t)) = 0} for all i,j 

Ai (?(<)) > for all i,j 

ff (0) = q° such that Aj(<?°) > for all i,j , q(0) = u° 
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The main equation 

Mq-f = J2 XijGijiq) G -N(Q , o) (45) 

i<j 

expresses the fact that overlapping is prevented by a repulsive force (the impulsion) acting on 
each sphere along the normal vector at the contact point. When there is no contact, N(Qq, q) is 
reduced to {0}, so that (|45p reads as Mq = f, which is the Fundamental Principle of Dynamics 
applied on each sphere. Equation q(t + ) = Pc q(t) <i(t~) provides the inelastic collision model. It 
can be extended to an elastic collision model with a restitution coefficient e by writing 

q(t + ) = P Cg(t) q(t~) - eP N(QoMt)) q(t-). 

We assume for simplicity that each mass m« is equal to 1. Then Problem (|44p fits into the 
previously studied framework. 

Remark 4.1 The case of different masses can be taken into account by using the adapted scalar 
product (u,v)m = (Mu,v), as done in JSj/. It turns back to replace the projection step in the 
numerical algorithm by 

U n+1 =V Kh{t ^, q n h) (u n + hM-lr), 

where P here denotes the projection relatively to this new norm. 

M being a diagonal matrix with non-negative diagonal coefficients, it is easy to show that the 
following results still hold true in that case. 

We emphasize that Assumption (|A0p is satisfied as soon as 

minr,- > 0, 



and then Assumptions (|Aip and (|A4p hold true. 

In order to apply our previous results, it remains to check Assumption (|A6p . As explained in 
[28] . that corresponds to estimate the Kuhn- Tucker multipliers. Such an estimate is given in the 
following lemma. 

Lemma 4.2 There exists a > (depending on N and on the radii ri) such that for all q G R 3Af , 
F G M 3JV and Lagrange multipliers (/%) G f^t^ -1 )/ 2 satisfying 



fiijGij(q) = F with /ly > and pij = when Dij(q) > 0, 



then 



Hij < a\F\. 

Concerning the proof of this lemma, we refer the reader to Proposition 4.7 of [28] (for a geometric 
proof) and to Proposition 2.18 of [U] (for a more "physical" proof). These proofs are written in 
a 2-dimensional framework but they can be easily extended in our 3-dimensional case. Actually, 
Lemma 14.21 is equivalent to Assumption (|A6p with p = 0. However, it can be extended and 
still holds for p small enough (for example p < infjrj), see Remark 4.11 of [28]. Consequently, 
Assumption (|A6p is satisfied for some small enough p > 0. 

According to our main Theorem (Theorem 1 1.3j) . it follows that Problem (|44p has solutions and 
that the associated numerical scheme converges (up to a subsequence). We can allow the radii 
to depend on time as soon as is uniformly twice-differentiable in time and 

inf inf rAt) > 0. 
te[o,T] i 

These theoretical results permit to legitimate the implementation of this numerical scheme. This 
was performed by the second author by creating SCoPI Software [27] . We refer the reader to [8] 
for some good properties of stability and robustness for the algorithm and efficiency for large 
time steps. 
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Remark 4.3 We refer the reader to [1], where the second author extends this model in order to 
consider gluey particles. In this case, she add an extra parameter (depending on q) for describing 
the corresponding admissible set. This new operation does not keep the necessary regularity of 
the admissible set. She has already obtained a result of convergence for the associated numerical 
scheme in the single- constraint case. We plan in a forthcoming work to extend this proof with 
the ideas presented here in order to deal with the multi- constraint case. 
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